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Abstract 


High-lift system aerodynamics has been gaining attention in recent 
years. In an effort to improve aircraft performance, comprehensive studies of 
multi-element airfoil systems are being undertaken in wind-tunnel and flight 
experiments. Recent developments in Computational Fluid Dynamics (CFD) 
offer a relatively inexpensive alternative for studying complex viscous flows by 
numerically solving the Navier-Stokes (N-S) equations. Current limitations in 
computer resources restrict practical high-lift N-S computations to two 
dimensions, but CFD predictions can yield tremendous insight into flow 
structure, interactions between airfoil elements, and effects of changes in airfoil 
geometry or free-stream conditions. These codes are very accurate when 
compared to strictly 2-D data provided by wind-tunnel testing, as will be shown 
here. Yet. additional challenges must be faced in the analysis of a production 
aircraft wing section, such as that of the NASA Langley Transport Systems 
Research Vehicle (TSRV). A primary issue is the sweep theory used to 
correlate 2-D predictions with 3-D flight results, accounting for sweep, taper, 
and finite wing effects. Other computational issues addressed here include the 
effects of surface roughness of the geometry, cove shape modeling, grid 
topology, and transition specification. The sensitivity of the flow to changing 
free-stream conditions is investigated. In addition, the effects of Gurney flaps 
on the aerodynamic characteristics of the airfoil system are predicted. 
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1 . Introduction 


The performance of modern transport aircraft, in terms of payload and 
range, is governed by the performance of the high-lift system. The increased lift 
coefficients required during low speed flight, associated with take-off and 
landing approach, are achieved by deploying leading-edge slats and trailing- 
edge slotted flap systems. Design objectives include improved lift and drag 
characteristics for a given weight and complexity, or decreased weight and 
complexity for a given Cl requirement. The combination of airfoil elements in a 
multi-element system can produce remarkably better performance than the sum 
of individual contributions, and this synergistic interaction involves complex flow 
physics, making analysis and design difficult. The lack of understanding in this 
area has inspired numerous comprehensive testing programs, both in wind 
tunnels and in actual flightT2. Unfortunately, the expense and time involved in 
instrumentation and data collection severely limits the rr»nge of test results. 
Furthermore, wind tunnel testing can rarely be accomplished at full-scale 
Reynolds numbers and the extrapolation to flight is non-linear due to various 
Reynolds number effects^ (Fig, 1.1). 

A potential solution to the challenges of testing programs emerged with 
the computer revolution of the 1980's and the development of Computational 
Fluid Dynamics (CFD) methods. Computer simulations of flows were cheaper 
and much faster than experimental programs, so more configurations could be 
analyzed, and full-scale flow parameters could be applied. The rapid growth of 
CFD methods and computer power lead to speculation that CFD would soon 
render the wind tunnel obsolete. However, ail computational methods suffer 
from limitations in terms of simplifications or assumptions made in the 



mathematical model, or hardware limits in CPU power, speed, and memory. 
Early CFD tools used potential flow and panel methods which solved for 
inviscid, irrotational flows. Later, Euler methods were applied, allowing regions 
of rotational flow to exist in the solution. Viscous effects were modeled by 
coupling inviscid solvers with boundary layer methods. 

Recently, several codes have been developed that solve the Navier- 
Stokes equations, the most accurate mathematical description of flow known to 
date. Yet, even these codes suffer from assumptions about the flow and 
computational weaknesses. Grids used to discretize the flow field are generally 
not fine enough to resolve all the details of a true flow. Current codes solve the 
Reynolds averaged form of the Navier-Stokes (N-S) equations, in which the 
effects of small scale turbulence are approximated by a turbulence model. 
Furthermore, artificial viscosity is required In order to achieve numerical 
convergence. Finally, due to computer limitations in both memory and CPU 
speed, most N*S high-lift studies are 2-D. Assuming 2-D flow, in turn, requires 
modeling assumptions to account for sweep, taper, and finite wing effects. 
Despite the limitations, these N-S methods provide the most complete viscous- 
flow analysis available. 

The purpose of this study is to analyze the existing multi-element wing of 
the NASA Langley Transport Systems Research Vehicle (TSRV) using 2-D 
Navier-Stokes methods and correlate with wind tunnel and flight test data. 
Following this introduction is a brief review of multi-element airfoil flow issues in 
Section 2. Section 3 provides a description of the computational tools used in 
this study, along with 2-D test-case results. Section 3 ends with a description of 
the sweep theories which are applied to correlate with 3-D flight-test data. 
Section 4 provides a description of the flight experiment of the TSRV, which 
provides the basis of this study. Section 5 steps through the effects of surface 



roughness, cove shape, sweep theory, grid refinement, and transition 
specification, leading up to the comparison with flight results. Next, the ability to 
predict the effects of changing free-stream conditions, in terms of angle of 
attack and Reync'ds number are addressed. Finally, an analysis of the addition 
of Gurney flaps to the aft flap are shown. Section 6 summarizes the 
conclusions of this research. Throughout this text, various suggestions ' ' 
future improvements will be made when appropriate. 
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Figure 1.1 • Scale effects on maximum lift (From Reference 1). 


2. Multi-Element Airfoil Flows 


The effects of slots in multi-element airfoils were highly misunderstood 
before A.M.O. Smith's classic papers on high-lift aerodynamics^-®. Previously, 
the commonly accepted explanation for slotted airfoil performance 
enhancement was that the slot allowed high energy lower-surface flow to re- 
energize the upper surface boundary layer, thus delaying flow separation. On 
the contrary, Smith pointed out that slot flow is low in energy and the 
advantages of slotted airfoils involve mostly inviscid phenomena. Achieving 
improved high-lift performance requires understanding the boundary layer and 
separation, and finding the inviscid pressure distribution which produces the 
least stress on the boundary layer, thus reducing the tendency for flow 
separation. 

As described by Smith, there are five predominant favorable effects of 
gaps in multi-element airfoil flows. The circulation of a forward element induces 
flow on a trailing element counter to the natural acceleration around the leading 
edge. This slat effect reduces the leading-edge suction peak on the trailing 
element, thus reducing pressure recovery demands and oelaying separation. 
The trailing element, however, induces a circulation effect on the forward 
element which tends to increase the loading on the forward element, increasing 
the lift, but also increasing pressure recovery demands. Yet, the high velocity 
flow on the upper surface of the trailing element allows the flow to leave the 
forward element at a higher speed. This dumping effect reduces the pressure 
recovery of the forward element and favors off the surface pressure 
recovery, which is known to be more efficient than recovery in contact with a 
wall. Finally, each element has a fresh boundary iayer which originates cn 



that element. A thin boundary layer can withstand stronger pressure gradients 
than a thick one and is less likely to separate. Effectively, the overall pressure 
recovery of the multi-element system is split by all the elements, but the 
boundary layer does not continuously grow along the chord as it would if the 
system was a single element. 

The primary viscous effect of slots is the existence of individual wakes 
from each element of the system. These wakes are thought to provide a 
damping effect on the pressure peak of trailing elements, reducing the tendency 
of the flow to separate^. Yet, the wakes often tend to merge with the boundary 
layer of the trailing element. The resulting confluent boundary layer is much 
thicker than an ordinary boundary layer, so the likelihood of separation 
increases. Clearly, gap size optimization requires a balance between the 
inviscid and viscous effects which favor smaller and larger gaps, respectively. 



3. Computational Methods 


Multi-element flow solutions entail complex inviscid and viscous flow 
phenomena due to gap and cove structure and wake interactions. However, 
many conventional computational studies, particularly those used in design, use 
inviscid solution procedures with limited or no boundary layer modeling. While 
providing general insight into the flow development, such methods neglect 
many issues of importance, such as confluent boundary layers, cove flows and 
upper surface separation, and transition. While 3-D viscous grid generators 
and solvers are becoming available, current computational limitations, both in 
memory and CPU time preclude practical application. A logical compromise is 
a detailed 2-D viscous solution which requires computational power of the same 
order as 3-D inviscid solutions. The current study employs two currently 
available 2-D unstructured grid Navier-Stokes solvers. 

3.1. Grid Generation 

The flow field surrounding a multi-element airfoil is discretized with an 
unstructured g^d consisting of triangular cells. This approach offers flexibility to 
conform to complex geometry and adapt to flow solutions. Cells are easily 
added where local refinement is needed to capture high gradients in the flow 
variables. Unfortunately, the lack of natural connectivity which gives 
unstructured grids such flexibility causes additional memory overhead. 
Alternative structured-grid approaches, such as the Chimera overset method® 
and the block structured approach^ benefit from implicit connectivity within each 
zone, but suffer at interpolation interfaces and experience grid "striping” if 
adapted. Presently, no method ic clearly superior end only the unstructured 
approach is discussed here. 



With the oDjective of capturing viscous phenomena in the flow solution, 
the importance of grid topology cannot be overstated. In order to resolve the 
high gradients involved in shear flows, dense grid cell distribution is necessary. 
Yet, global refinement results in wasted memory and computation time in the 
areas where not needed. Thus, "efficient" grids utilize stretching factors to vary 
cell size between viscous and inviscid regions of the flow. Furthermore, viscous 
gradients commonly require much higher resolution in the normal direction than 
in the tangential direction. Figure 3.1 illustrates how high aspect ratio grid cells 
can be used to resolve viscous layers without an excessive number of nodes. 
While these high aspect ratio cells degrade the local efficiency and accuracy of 
the flow solver, the increase in global efficiency justifies the use of such cells. 

Two methods of unstructured grid generation were applied. The first is 
part of the NSU2D® package developed by Dimitri J. Mavriplis under the support 
of ICASE and McDonnell Douglas Aerospace. This method is designed to 
provide high aspect ratio viscous spacing near solid surfaces as well as in the 
wake regions in order to resolve boundary layer and wake flow structure. The 
other grid generation tool considered here is USVGRID^ writ‘en by Shahyar 
Pirzadeh under the support of the NASA Langley Research Jenter. This 
method provides viscous type cells at solid surfaces and distributes cells 
smoothly around even the most complex shapes. 

3.1.1. NSU2D package 

Grid generation in the NSU2D package involves two primary steps: 
node generation, then tnangulation. The point distribution comes from an 
interactive module called Gridpts. These nodes span the domain to be gridded, 
with coarse spacing in the outer region and fine spacing in the vicinity of airfoil 
surfaces. The output of Gridpts is triangulated in batch mode by UMESH2D. 
The topology of the resulting grid is largely determined by the point distribution 



used for the iriangulation. 
3.1.1 .1. GridPts 


The Gridpts module provides an interactive interface to distribute node 
points for the grid. The multi-element airfoil geometry is input. The user selects 
break points between which piecewise splines are fit to the geometry. A quick 
panel solution of flow past the configuration determines the approximate 
location of the wakes. The splined body shape and wake of each element 
forms an inner boundary for a structured C-mesh around that individual 
element. Normal spacing at the body, the number of points around the body, 
and the distance to the outer edge of the grid are among the inputs to the 
structured hyperbolic grid generator. The outer boundary of the main element 
structured grid will be the outer boundary of the final unstructured grid around 
the multi-element system. Note that up to this point, the process could be used 
for a Chimera overset grid method, but now we deviate. 

The set of all nodes contained in the structured grids must be 
manipulated to provide a suitable point distribution for a single unstructured grid 
around the system. The structured grids are coarsened as desired to be used 
for multigrid convergence acceleration in the flow solver. Any points that fall 
within the body of any element are eliminated. Interactive sessions allow the 
specification of a "preserved region" around each element where the structure 
of the grid nodes is to be exactly preserved in the final triangulated grid. The 
final operations of Gridpts are a smoothing of the node distribution and removal 
of excess points. In particular, specification of a maximum cell aspect ratio and 
distance downstream to achieve it dictates the blending of viscous wake ceils 
into the downstream region. Similarly, a maximum aspect ratio in the direction 
normal to the elements smoothes the inviscid point distribution. 

Gridpts writes many intermediate files during the procedure. Thus, the 



process may be interrupted and resun.ed at these intermediate points. Log files 
further expedite operation by providing user responses to various input needs 
(while preserving the interactivity of the graphic sessions). Log files are 
particularly helpful for the creation of multiple grids of varying coarseness, that 
are otherwise very similar, to be used with a multigrid algorithm or for a mesh 
refinement study. The final output file of Gridpts contains all the information 
needed in UMESH2D. It includes the set of nodes, the extent of the preserved 
regions, and stretching factors based on the aspect ratio distribution of the 
structured meshes. 

3.1.1.2. UMESH2D 

The points to be triangulated consist of two regions. The "inner region" 
incorporates the set of preserved structured grid regions around the elements. 
This region is triangulated simply by connecting a diagonal across each four- 
sided cell, so the original spacing is preserved. The "outer region" 
encompasses the rest of the domain and requires special treatment in order to 
achieve a smooth variation of cells from the high aspect ratio viscous region to 
the low aspect ratio inviscid region. The entire process is carried out in batch 
mode ^uided by a small input file. 

Triangulation of the outer region is accomplished using a Delaunay point 
insertion process’^. The basic algorithm, credited to Bowyer^^ provides a 
method to systematically insert the desired node points into the triangulated 
region, ‘^he initial triangulation consists of cells extending from the inner to 
outei boundaries. As each point is inserted, all neighboring triangles with 
circumcircles that contain the new point are flagged. The edges of these 
existing triangles are removed and the resultant cavity region is re-triangulated 
by connecting the new node to every node along the edge of this cavity. The 
resulting mesh is a Delaunay triangulation of the given point set. This 



construction has certain well-known properties, including a tendency to contain 
low aspect ratio cells. While acceptable for the inviscid region, this outer grid 
will not blend well with the high aspect ratio cells in the inner region, and abrupt 
changes in cell shape can be damaging to the solution. Thus, edge swapping 
operations revise the triangulation using the stretching vectors taken from the 
origina; structured meshes. The nodes are mapped to a stretched space for 
triangulation. When mapped back to physical space, the cells will flatten in 
proportion to the stretching factor, resulting in a smooth distribution of aspect 
ratio through the grid. 

The resulting grid may be viewed using MESHPLT, which is available 
with the NSU2D package. Careful examination of the grid topology reduces the 
risk of failure in the solution process. 

3. 1.1. 3. Limitations 

While providing the means for generating excellent grids for viscous 
calculations, various limitations are addressed. As perhaps inferred in the 
preceding description, the process of generating grids using Gridpts and 
UMESH2D can be long and tedious, requiring several iterations before the 
desired topology is achieved (though log files partially relieve this effort). 
Typically, each iteration takes about fifteen minutes to perform, so generating a 
satisfactory grid may take on the order of a couple of hours. Also, the inherent 
ability of unstructured grids to fit any geometry is limited by the intermediate use 
of structured meshes in the process. Robustness problems arise in the 
structured grid generation, particularly in concave regions such as coves. To 
address this problem, the option to use alternative structured mesh generators 
is included in Gridpts. The value of the structured meshes in the process 
justifies the difficulties in their use. In addition to robustness problems in 
Gridpts, UMESH2D suffers from round-off errors in the point Insertion process. 



though these can usually be corrected by varying UMESH2D inputs. Another 
limitation of this process is the difficulty in achieving a smoothly blended mesh 
when long regions of high aspect ratio ceils are retained in the wakes of the 
high-lift elements. Decreasing the extent of the viscous wake spacing 
smoothes the mesh but causes additional dissipation in the solution which 
reduces the ability to capture certain viscous interactions such as confluent 
boundary layers. Furthermore, generating smooth grids through thin gaps is 
difficult due to the preservation of structured mesh regions with high aspect 
ratio cells. Generally, as with any process, trade-offs must be made, and 
problems in grid generation with the NSU2D package can be averted with 
experience. 

3.1.2. USVGRID 

USVGRID provides a streamlined approach to unstructured grid 
generation. A single input file contains all the required information for the grid 
generation. USVGRID operates non-interactively with an optional graphical 
interface to view intermediate steps in the process, or can be run as a batch 
job. Euler grids, or Navier-Stokes grids with viscous spacing near the airfoil 
surfaces may be produced. Two common output formats are available so the 
grid may be used with a variety of available solution codes. 

Qp.er.atiOQ 

The primary hurdle in using USVGRID is preparation of the input file. A 
numbered list of points defines the surfaces of the airfoil geometry and outer 
boundaries. Piecewise splines through the given geometry points are indicated 
in lists of all the points co nposing each spline. Next is some information about 
the background mesh, followed by the locations and strengths of 'sources* 
which govern the size and clustering of the ceils in the mesh. The region of 
influence for each source is also specif d. Finally, for viscous grids, the normal 



spacing at the body and a stretching factor to control the rate of decrease of cell 
aspect ratio in the normal direction are listed. 

The generation of viscous grids with USVGRID is accomplished with an 
Advancing Layer/Advancing Front approach^. The grid cells are created as the 
front marches outward from the inner boundaries. The viscous grid for the 
region near solid surfaces is generated by the Advancing Layer technique in 
which nodes are added along lines normal to the surface using the specified 
spacing and stretching. Beyond a certain distance from the surface, the 
Advancing Front technique is used to provide an isotropic triangulation of the 
remaining inviscid region. Grid "shocks"^ 2 ^ which typically occur in structured 
grids near concave corners, are avoided in this unstructured approach by a test 
procedure which predicts when cells in the current march will overlap and 
reduces the number of new faces on the next front accordingly. The entire grid 
generation process is rapid, one iteration taking about five minutes. The 
resulting grid is very smooth, generally with more gradual variation in cell size 
through the mesh than in grids generated with UMESH2D, 

3. 1.2.2. Limitations 

The primary difficulties involved in USVGRID are associated with the 
complex input file. Specification of the spline lists and source terms will be 
much simpler when a graphical interface for interactive specification becomes 
available; such an interface is under development. Until available, use of a 
pre-processing code is prudent. Placement and properties of source terms 
require several levels of iteration until the desired node distribution, particularly 
near the solid boundaries, is achieved. Then, when the desired resolution is 
obtained near the airfoils, the grid tends to have excess nodes away from the 
airfoils and will require more computation time in the solver. Also, no high 
aspect ratio cells are generated downstream of the elements for fine resolution 



of wakes and confluence. Perhaps an option could be to define a wake line 
along with the geometry from which to march cells Initially using the advancing 
layers approach used on the element surfaces. 

3.2. Mesh Adaptive Refinement 

Considering the influence of grid topology on the flow solution, the ability 
to refine existing grids to better suit a particular solution is considered. As 
mentioned previously, local refinement to capture high gradients is more 
efficient than global refinement. The unstructured grid generation process 
distributes cells with refinement where the need is expected. However, without 
a priori knowledge of the solution, the grid will likely lack the needed resolution 
in some areas of the flow. For instance, wakes will dissipate where the grid is 
too coarse, so the occurrence and effects of confluent boundary layers cannot 
be addressed. Furthermore, wake locations change in different operating 
conditions, so grid clustering may not exist where needed. Also, failure to 
accurately predict pressure peaks can lead to entropy generation and other 
globally damaging effects in the solution. The motivation for solution-adaptive 
refinement is the ability to selectively add nodes only in regions where gradients 
in the flow variables are large. The refined mesh either replaces the existing 
mesh for a new calculation or may be added to the top of a multi-grid sequence. 

3.2.1. NSU2D Adaptive Refinement 

The NSU2D package includes MESHAD for solution-adapted grid 
refinement. MESHAD uses the geometry spline information (an intermediate 
Gridpts file), the existing finest grid, and a restart file from the NSU2D solver to 
generate a new grid. In the MESHAD input file, each flow variable to be 
examined and the criterion for adding nodes is specified. The flow variables 
available for adaptation are density, pressure, velocity magnitude, and Mach 
number. Resolution of inviscid phenomena, such as pressure peaks can be 



achieved using density or pressure, whereas high Mach number and velocity 
gradients occur in regions of viscous shear flow. Any combination of these 
variables may be used. The criterion for refinement to a particular property is 
the ratio of the local change in that property within one cell to the average 
change within all cells in the domain. The allowable ratio is specified in the 
input file, and any cell that exceeds this ratio has new points added at the 
midpoints between the existing nodes. Specifying the ratio as 0 results in 
global refinement, producing four times as many cells in the grid, and a ratio of 
1 results in refinement anywhere the gradient is higher than average. Once all 
new points are determined, they are added to the existing mesh using the point 
insertion process as in UMESH2D. 

3.2.2. USVGRID 

No solution-adapted refinement code is currently available for use with 
USVGRID (partially attributed to the fact that development of USVGRID is 
independent of a particular solution code). Manual refinement can be 
accomplished by viewing a solution and adjusting the source placement and 
strengths in the input file and generating a new grid better suited to the solution. 

3.3. Navier-Stokes Solvers 

The full Navier-Stokes equations provide a complete mathematical 
model for flow solutions. However, computational limitations in both hardware 
and software require assumptions to simplify the process. In order to reduce 
computer processing and memory requirements, 2-D flow is studied here. Still, 
the length scale of turbulent phenomena is generally smaller than the cell size 
of current grids, so the equations are Reynolds averaged and the mean effects 
of the small scale turbulence are accounted for by a turbulence model. The 
purity of the mathematical model is further compromised by the requirement of 
numerical stability, which is achieved by adding artificial viscosity to the 



equations^ 3. Available solvers are utilized in this study, but the understanding 
of the solution is incomplete without knowledge of the assumptions and 
limitations in the formulation. 

3.3.1. Governing Equations 

The Reynolds-averaged Navier-Stokes Equations govern the continuity 
of mass, momentum, and energy, and are given as; 

(3-1) 

where Re» is the free-stream Reynolds number based on the free-stream 
velocity, U« and the reference chord, c, and w is the set of conserved variables: 



with density, p, velocity, 0 = ui +v] , and total energy, E. The convective 
fluxes, F^, consist of algebraic functions of the conserved variables and the 

pressure, which is related to the conserved variables through the equation of 
state for a perfect gas. The viscous fluxes, F„, consist of functions of the first 

derivatives of the velocity components, as well as molecular viscosity, p, and 
turbulent eddy viscosity, pt- (For details, consult References 8, 19). 
Sutherland's Law may be used to compute p, but pt must be solved for 
simultaneously with the conservation relations using the additional equation(s) 
provided by a turbulence model. In this formulation, the global effects of small 
scale turbulence are imposed on the solution without requiring the extreme grid 
resolution which would otherwise be necessary to capture the full turbulence 
effects. Consequently, the quality of the solution depends on the ability of the 
turbulence model to accurately predict turbulence levels which contribute to the 



viscous stress and heat «lux in the flow. The turbulence model used in this 
study is the one-equation model of Spalart and Allmaras’^, which has been 
found to provide better resuiis than many alternative models^s-’®. 

3.3.2. NSU2D 

The steady-state flow solution is calculated by NSU2D using an explicit 
Runge-Kutta time stepping scheme applied to a finite element discretizations. 
Convective terms are central differenced and two forms of artificial viscosity are 
used: second-order accurate biharmonic operators throughout the domain, and 
first-order accurate Laplacian terms in regions of high gradients or shocks. As 
with any explicit scheme, the time step used is limited by a stability condition 
known as the CFL number, named after Courant, Friedrichs, and Lewyis The 
exact form of the CFL number depends on the discretization, but is generally of 
the form: CFL = |c^|, with CFL typically less than 1 for stability of a single-stage 

explicit scheme. The five-stage Runge-Kutta integration scheme applied in 
NSU2D allows the use of CFL numbers up to around 3.5^T 

Since most studies concentrate on the steady state solution and not the 
transient, the definition of "time" is relaxed, allowing various methods of 
convergence acceleration to be used. A simple method is to observe the 
dependence of a stable time step on the cell size and thus allow the time step 
to vary throughout the mesh. In addition to local time stepping, residual 
smoothing can further accelerate convergence. By implicitly averaging 
residuals with those of neighboring nodes, the stability limit of the scheme 
increases, allowing the use of larger CFL numbers. A CFL number of 5 was 
found to provide a good balance between efficiency and robustness. 

Another more recently developed and more powerful convergence 
acceleration technique is multi-griddingi®. Solutions rapidly obtained on 
sequentially coarser grids provide improved corrections for the fine grid, on 



which detail in the flow solution is obtained. Typical multigrid cycles are 
illustrated in Figure 3.2. The larger cells oi H'e coarse meshes allow more rapid 
travel of boundary information, as well at *■ i^voing out high frequency noise in 
*' 'olution. Independent triangular metht., lake up the multigrid sequence 
a her than using subset coarsened ■ or polyhedral agglomeration^ 9 
techniques, allowing control over giid . . i . gy and allowing arbitrarily refined 
mesnes. The transfer of inform^; oc t ..een meshes is accomplished by a 
simple linear -n.erpolation br hwr . lode of one mesh and the nodes of the 
surrounding cel! on the other m ; 3h. This information is stored as 3 interpolation 
addresses and 3 weights for <rv-ich node, as determined in a pre-processing 
operation. 

Operation of NSU2D involves a primary input file which provides details 
for the solution process, along with the names of the grid files and if any laminar 
flow regions are to be specified, a boundary-layer transition file. The flow solver 
may be compiled and run on a supercomputer or a workstation if sufficient 
memory is available. Typically between 10 to 20 Mwords of RAM on a Cray, or 
from 64 to 96 Mbytes of RAM on a workstation are needed for the 
computations, depending on the grid size. Upon completion, a restart file 
containing all the flow variables at the current time step is written, allowing the 
job to be re-submiUed if further convergence is desired. The solution for the 
entire flow field is written to a file which may be viewed using PLOT2D, whicii is 
available with NSU2D. The finest grid, contours of various flow properties, and 
velocity vectors may be viewed. A post-processing operation may be used to 
extract velocity profiles comparable to bou*^dary layer rakes used in 
experimental studies. Another file contains convergence information and 
surface values which may be post-processed to plot surface properties such as 
pressure coefficiet.' 



3.3.3. FUN2D 

FUN2D was developed at NASA/Langley Research Center by W. Kyle 
Anderson20. An implicit, upwind finite element scheme is used to march the 
solution to steady state. The implicit formulation provides numerical stability, 
allowing higher CFL numbers, and thus larger time steps, to be used. FUN2D 
is stable using a CFL number on the order of 100, as compared to around 5 for 
NSU2D. The apparent gain in speed is balanced b', the longer time per 
iteration than the explicit scheme due to the need for solving a large 
simultaneous system of algebraic equations at each time level, resulting in 
comparable times for numerical convergence between the codes. Artificial 
viscosity is not explicitly added in the equations, out the stabilizing effect is 
provided by the upwinding of the convective terms in the equations (See 
Appendix 1). Convergence acceleration techniques such as multigridding and 
local time stepping are implemented as in NSU2D. 

Operation of FUN2D is much like that of NSU2D, including a restart 
option and comparable output files. The output is formatted for use with the 
FLANAL visualization package developed by Kelvin Edwards at NASA/Langley 
Research Center. As will be shown in the lollowing section, the results of 
NSU2D and F' iN2D are in good agreement for a three-element configuration, 
as they are for all configurations studied. Thus, only a brief comparison 
between solvers is conducted, and most solutions are obtained using NSU2D. 

3.4. 2-0 Test Cases 

Before discussing 3-D flow fields, an evaluation of the numerical 
methods using strictly 2-D test cases is advantageous. Computational results 
include error introduced by the formulation, discretization, and computer round- 
off. In addition, the correlation with 3-0 data involves error due to sweep 
corrections, as well as the lack of spanwise flow, interference, and wing-tip 



e.fects in Ihe 2-D calculations. Furthermore, actual aircraft used in flight exhibit 
more of a tendency to deform aeroelastically than wind-tunnel models. Thus, 
comparison with 2-D wind-tunnel data greatly reduces the sources of potential 
error and allows an objective evaluation of the capability of the codes. 

The Low Turbulence Pressure Tunnel (LTPT) at NASA Langley provides 
an abundance of reliable test data including several tests of multi-element wing 
sections at Reynolds numbers appi caching those of flight. The first case 
considered here is a three-element airfoil available from McDonnell Douglas 
Aerospace which was used as a benchmark test case in the High-Litt Workshop 
at NASA Langley in May 19932’. The second is a four-element airfoil 
developed by the Boeing Commercial Airplane Group. The geometry of this 
configuration remains proprietary. 

3.4.1. Three Element 

The three-element airfoil tested is McDonnell Douglas model MDA 
LB-546, which consists of a slat, main element, and a flap. The particular case 
used for this comparison is Geometry A of the High Lift Workshop, which is a 
landing type configuration with a 30° flap deflection. Free-stream conditions are 
characterized by a Mach number of 0.20. Reynolds number of 5 million, and 
angle of attack of 8.12°. A typical grid used for the calculation is shown in 
Figure 3.3. This grid contains 48121 nodes, with 291 nodes on the main 
element surface, 319 nodes on the slat, and 265 nodes on the flap. Two 
coarser grids, each with about one fourth as many nodes as the previous, are 
used for the multi-grid sequence. The NSU2D solution requires about 8.5 
Mwords of memory and 30 minutes of CPU time on a single processor of a 
Cray Y-MP supercomputer. FUN2D requires about 28 Mwords of memory and 
25 minutes of single-processor Cray Y-MP time. The required memory per 
node is approximately 180 words/node for NSU2D and about 600 words/node 



for FUN2D. On an IRIS Indigo 2 workstation, using double precision, an 
NSU2D solution requires around 64 Mbytes of memory and on the order of 10 
hours of CPU time. The convergence histories for NSU2D and FUN2D are 
shown in Figures 3.4 and 3.5, respectively. Both plots show the decrease of a 
local quantity, namely the maximum rate of change of density or mass, and the 
asymptotic convergence of a global quantity, the lift coefficient. For both codes, 
a four order of magnitude reduction in the residual, and a convergence of the lift 
coefficient to the fourth decimal place is accomplished within one hour of Cray 
CPU time for all grids tested. 

The surface pressure distribution is well predicted by both codes, which 
are in close agreement, as shown in Figure 3.6. The overprediction of the slat 
pressure may be due to poor representation of the separated cove flow due to 
limitations of the turbulence model, resulting in increased circulation on the 
element. The overprediction of the main element pressure peak may also 
contribute to the ovet prediction of the slat loading through the circulation effect. 
The flap pressure comparison is very good. As previously mentioned, since 
NSU2D and FUN2D give very close to the same results, no further comparisons 
between the two codes are shown. 

While the surface pressure comparison shown here resembles the state 
of the art seen at the workshop and in the literature, the ability to better capture 
the flow structure is addressed with a grid resolution study. The grid presented 
above is used as the baseline case. A coarse grid comparison simply uses a 
solution obtained on the next grid of the multigrid sequence used in the baseline 
case. A finer mesh was also generated to test grid dependence, and was also 
designed to include long regions of high aspect ratio cells in the wake region of 
each element in orde.* to discern the value of including such cells. The sizes of 
the three meshes are listed in Table 3.1. 
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Baseline Mesh 

Coarse mesh 

Fine Mesh 

Total Nodes 

48,121 

12,173 

70,720 

Slat Surface 

319 

160 

299 

Main Element Surface 

291 

147 

379 

Flap Surface 

265 

133 

299 


Table 3.1 - Mesh sizes for three-element grid comparisons. 


The computed surface pressure distributions for the three mesh cases 
are shown in Figure 3.7. Perhaps even more revealing are boundary-layer 
profiles, which are shown in Figure 3.8. The experimental profiles were taken 
using a boundary-layer traverser. The computational results are found by 
computing a surface normal at each chordwise station desired and taking the 
component of velocity perpendicular to this line by interpolating between nodes 
of the grid. The coarse mesh clearly ha? normal spacing near the surface 
which is too laige to resolve the boundary layers. The artificially thick boundary 
layers and wakes tend to dampen the pressures on trailing elements and thus 
reduce the circulation of the entire system. The pressure prediction of the 
baseline and fine mesh are very similar, suggesting that ihe baseline mesh is 
adequate for general solutions and the additional high aspect ratio cells carried 
into the wake are not necessary for prediction of the pressures. 

The velocity profiles are also very similar for the baseline- and fine-mesh 
solutions, both resolving wakes much better than the coarse-mesh solution. 
The extended region of high aspect ratio cells in the tine grid does reduce wake 
thickness slightly and improves the velocity levels. The fine grid solution shows 
a slat wake at the r\ = 0.45 station. The experimental result has only a slight 
perturbation at this location, but the experi.mental slat wake may have been 




disturbed by the instrumentation and tunnel turbulence, and should be 
somewhat thinner than the fully-turbulent computations which do not include 
any laminar flow on the slat. Dissipation in the computational solutions due to 
cell size and artificial viscosity is evidenced by the spreading of the wake as it 
convects downstream. With the coarse grid, the slat- and main-element wakes 
completely dissapear within a short distance downstream. Near the trailing 
edge of the flap (r\ = 1.122;, there is a discrepancy between computation and 
experiment on the location of the wake center. The incipient separation in the 
experimental profile is not predicted in the computations, suggesting the need 
for further development of turbulence models. The multi-element flow structure 
is best resolved with fine grid spacing in the wake regions, so the calculations in 
this study are made using grids generated with the NSU2D package rather than 
USVGRID. 

3.4.2. Four Element 

While the three-element results show excellent agreement between 
computations and experimental data, several issues remain to be addressed 
before comparison with the five-element results of the TSRV. First, the 
accuracy of the codes rrtay degrade as the flow field becomes more 
complicated by the addition of lifting-surface elements. Also, the interaction of 
multiple flap elements in close proximity is of interest. In particular, the four- 
element airfoil analyzed here includes a small vane element followed by the 
main flap element, not uniike the first two flap elements of the TSRV airfoil 
section. Because this geometry is a current interest of the Boeing Commercial 
Airplane Group, the geometry and results are not being published, so only a 
discussion of the results are included here. 

The solution was run for a flap setting appropriate for landing, at 5.9® 
angle of attack, a Mach number of 0.18, and three Reynolds numbers. 2.45-, 



6.70-, and 15.97-million. The finest grd for the solution contains 71812 nodes, 
with 253 nodes on the slat element surface, 371 nodes on the main element, 
301 nodes on the vane, and 299 nodes on the flap. The computational results 
from NSU2D and FUN2D closely match those of the wind-tunnel experimen . 
The computations also accurately predict the variation with Reynolds number. 
Both codes predict the rise in the leading-edge pressure peaks with Reynolds 
number, and the levels are accurate. There is a slight discrepancy in the slope 
of the pressure recovery near the flan trailing edge for the two higher Reynolds 
numbers. This discrepancy may be attributed to inaccuracies in the turb,.!v.;rice 
modeling in regions with high Reynolds stress, although no upper surface 
separation was evident. Overall, however, the results of this test case illustrate 
the ability of these codes to accurately predict the flow field and Reynolds- 
number effects for complicated 2-D multi-element configurations. 

3.5. 2-D to 3-D Correlation 

In order to apply 2-D methods to the study of aircraft wings, an 
understanding of the influences of 3-D effects on the flow is necessary. This 
study involves a 2-D analysis of an existing multi-element wing and correlation 
with flight-test data. In order to compare with the flight results, the effects of 
downwash, sweep, and taper on the 2-D solution must be accounted for. 
Currently, no method provides all 3-D factors short of a full 3-D solution, which 
is currently infeasible due to computer limitations (as discussed previously). 

The most basic method is known as a simple-sweep corroction^^. This is 
easily accomplished in pre- and post-processing operations, so modification of 
existing 2-D codes is not necessary. Simple-sweep theory is purely geometric, 
as illustrated in Fig 3.9. Applying simple sweep theory entails choosing a 2-D 
section of the wing in the freestream direction end rotating the cut by some 
nominal angle of sweep so the section is in a normal chordwise direction. 



Then, the 2*D solution is found and the results are corrected back to the tree- 
stream direction. Only the effect of sweep is accounted for, so this correction is 
only rigorously valid for an infinite-span, untapered, swept wing. 

In the pre-processing operation, the airfoil shape is modified by 
preserving a non-dimensional chord length so x/c is unchanged, and altering 
the y-values of the surface coordinates as: 

X = = (3_3) 

Cn VcJUj 

resulting in a thicker airfoil section (In reality, the thickness is unchanged and 
the chord is smaller). The components of the free-stream Mach number and 
Reynolds number in the normal direction are also found from the geometry: 

At = A/- cos A, (3-4) 


n (l>cosA)(^cosA) _ 2 . 

Re, = —2-2. = .i L = Re^cos^ A. 


(3»5) 

V V 

The 2-D solution is found for the modified section using Mn and Rep. Typical 
qualitative outputs such as the velocity vectors and pressure contours around 
the system are viewed in un-corrected form. Yet, any features that are to be 
directly compared with flight results must be converted back to the streamwise 
direction for comparison. Namely the surface Cp distribution, and any values of 
forces, such as the normal force, lift, and drag which were calculated before the 
sweep correction, must be "unswept". The procedure for correcting Cp begins 
with the calculation from the 2-D solution at each chordwise location, i: 


P - P 

Cp = — ! — a. 
' iYP.M, 


(3-6) 


Yet, we want Cp relative to the free-stream direction: 
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{yP^Ml 


(3-7) 


and this is achieved by multiplying Eq. (3-6) by the square of the ratio of the 
normal Mach number to the free-stream Mach number: 


Ca = 


P. - P- 


\yP^Ml 
and applying Eq. (3-4), this yields: 


(mA 

1 

(m ] 


= ^Pn, 


[mJ 

l^J 


(3-8) 


CPi=Cp„ (3-9) 

In addition to the simple-sweep rule, the correlation with flight data 
requires an angle of attack correction to account for the 3-D finite wing effect. 
Because the true wing has a finite span, wing-tip vortices and the associated 
downwash alter the direction of the effective incoming flow past the section. 
Unfortunately, the induced angle of attack is not known in general for the wing 
section, so the choice of an used in the 2-D computation is somewhat arbitrary 
and requires iteration. One ' .hod commonly used in conjunction with 
experiment is to choose an such that the section Cn matches that of 
experiment. However, this criterion is believed to be excessively optimistic, and 
the Cp distribution will certainty be close since the pressure integration yields 
the same normal force coefficient, Cn. The method of finding an used in this 
study is to choose that for which the slat element upper surface pressure 
distribution most closely matches the experimental level. The lead element is 
least likely to be affected by viscous effects due to its fresh, thin boundary layer 
and position upstream of the wakes of other elements. Using this dn> nwash 
correction, the computational inaccuracies in the downstream flow and in the 
prediction of forces on the system can mors realistically be arsessed. 

While simple sweep theory is easily implemented and provides relatively 
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good comparison with flight^a, it is subject to several limitations. One 
improvement to the simple-sweep theory is accomplished by accounting for 
wing taper. Rather than using a single nominal sweep angle for the 2-D to 3-D 
conversion, a local sweep angle is calculated for each chordwise station, i. 
Thus, the new y-values are calculated from; 


7 



JiL, 

cosA, ’ 


(3-11) 


Note that the only difference with the simple-sweep correction in Eq. (3-3) is 
that here the sweep angle, Ai, varies along the chord. The chordwise variation 
in sweep angle is linear in the tangent, as derived in Appendix 2. Thus, the 2-D 
geometry comes from a circular arc cut through the wing rather than a straight 
line, as sketched in Figure 3.10. Similarly, the Cp correction uses the local 
sweep, Aj. 


Based on the theory and knowledge of the wing being studied, some 
predictions about the local sweep rule results may be made. First, because the 
sweep angles are more accurate, the flow prediction should be more accurate. 
However, flap loadings tend to be overpredicted due to "blow-back" in flight and 
other geometrically-related errors, as well as grid resolution and computational 
errors. (Blow-back refers to the movement of the airfoil elements when loaded 
due to "play" in the flap linkages, reducing flap incidence angles and gaps.) 
Because the wing taper ratio is less than unity, and the wing sweep is rearward, 
the flap sweep angles are less than the nominal quarter chord sweep used in 
the simple sweep approach. Thus, from Eq. (3-9), we see that the error in tne 
simple sweep method will result in decreased Cp values over the flap elements 
compared to the local sweep approach, making simple sweep appear to 
correlate better with flight for the flap elements. 

A third method of sweep correction is based on the empirical results of 



the British "National High Lift Programme" reported by D. S. Woodward and D. 
E. Lean in Reference 24, and is implemented here as an addition to the local 
sweep procedure. In Reference 24, comparisons between wind tunnel results 
obtained with 2-D airfoils and wind tunnel results obtained with swept semispan 
wings revealed a consistent alteration to the Cp correction, of the form: 

Ca=Q„cos"A, (3-12) 

Rather than applying the exponent, n = 2, which is purely geometric, n was 
allowed to vary in the chordwise direction. The best correlation between 2-D 
and 3-D test data was found to use the exponent, n = 1.1 on the slat element, n 
= 2.1 on the main clement, and n = 3.6 on the flaps. The explanation for this 
variation in the value of the exponent is unknown. Yet the fact that the 
exponent n is near 2 for the main element agrees with the geometric theory. 
Perhaps the modification for the slat and flap elements partially accounts for the 
chordwise variation in downwash from the tip vortices of the 3-D model, which 
was observed in Reference 24. The higher downwash in the vicinity of flap 
elements decreases the effective angles of attack, thus reducing the loading. 
Since cosA < 1, the higher value of the exponent, n, for the flap elements will 
decrease the computed pressures of the 2-D solution, which will better match 
the 3-D pressures. Likewise, the smaller value of the exponent, n, for the slat 
will reduce the extent of the decrease of the pressure on the slat due to the 
sweep correction and will thus increase the an for the correlation. 

It should be mentioned that application of a different sweep correciion 
requires more than simply altering the numbers in pre- and post-processing 
operations. Since the geometry depends on the sweep rule, the entire 
computation process, including the generation of new grids, must be repeated 
for each new sweep correction. 
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a. Coarse grid; 9 nodes. 
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b. Globally refined grid; 81 nodes. 
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c. Stretched grid; 63 nodes. d. High aspect ratio grid; 21 nodes. 

Figure 3.1 - Grid topology at solid surface and typical velocity profile at /. 
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Figure 3.2 - Typical muitigrid cycles. 



Figure 3.3 <• Typical grid for three-element airfoil computations. 
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Figure 3.4 - Typical convergence history for an NSU2D solution. 



Figure 3.5 - Typical convergence history for a FUN2D solution. 




Figure 3.6 - Comparison between experimental and computed pressure 
distributions for a three-element airfoil. 

(M = 0.20, Re = 5 million, a = 8.12°) 



Figure 3.7 - Effect of grid topology on predicted pressure distribution 
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Figure 3.8 - Effect of grid topology on boundary*layer velocity profiles. 





4. Description of Flight Experiment 


The Transport Systems Research Vehicle (TSRV) is the prototype B737- 
100 which has been hiahly modified for high-iift flight research. This aircraft, 
shown in Figure 4.1, provides a good representation of typical transport 
systems in terms of size and planform, while the five-element airfoil provides a 
challenge for analysis methods and CFD code validation. Although future 
designs will likely involve less geometric complexity, the basic elements of 
multi-element flow physics are well represented by this test case. Test results 
indicate the existence of flow separation, attachment-line transition, boundary- 
layer relaminarization, and likely confluent boundary layers. 

The Phase I and II flight tests included a limited array of instrumentation 
to provide some preliminary insight into the flow, and key results were 
presented in Reference 1 . Pressure belts around each element at the 58% 
semispan station provided the chordwise load distribution. When flaps are 
deployed, these belts are not all at the same spanwise station, as shown in 
Figure 4.2, but the effect on the results is thought to be minor. Preston tubes 
placed on the main element upper- and lower- surfaces and the slat upper 
surface provided insight into the skin friction on the wing. The test matrix 
included a range of angle of attack for three pressure altitudes, five-, ten- and 
twenty-thousand feet. 

The upcoming Phase III flight tests scheduled for Spring, 1994 will 
include a more complete array of instrumentation to help answer many 
questions which have arisen in both the Phase I and II analysis and the current 
computational effort. Flush pressure ports will address the effect of pressure 
belts on the measurements. The pressure will be measured at five spanwise 



stations to test the two-dimensionality of the flow near the half semispan. A 
boundary-layer rake on each element will show boundary layer thickness, 
separation, and confluence. The Northern Digital OPTOTRAK^s system will 
reveal in-flight deflections due to blowback and aeroelastic deformations, 
helping to reduce geometry-related error in the analysis. Hot films will provide 
more detailed illustration of the boundary layer state. Finally, Gurney flaps of 
1% and 2% lengths will be tested. Decisions on the placement and sizing of the 
instrumentation were based on Phase I and II test results as well as 
preliminary computations. 
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Figure 4.1 - NASA Langiey TSRV test aircraft (B737-100). 
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Figure 4.2 - Phase II flight test Instrumentation iayout. 








5. TSRV Results 


The computational analysis of the TSRV wing section provides a means 
for evaluation of the capability of the current codes and sv/eep theories. The 
use of an actual aircraft wing and correlation with flight data reveal numerous 
difficulties encountered in a design environment which do not arise when 
comparing to 2*D wind tunnel data. In the following sections, various issues in 
the process leading up to the correlation with flight test results are addressed. 

A typical primary mesh for the five-element computations is shown in 
Figure 5.1. The spacing on the element surfaces and the distribution of nodes 
in the flow field is based on the grid properties used for the 2-D test cases , and 
is comparable to grids described in the literature. This mesh contains 72,975 
nodes, with 223 nodes on the slat, 399 on the main element, and 210, 313, and 
213 on the fore-, mid-, and aft- flap, respectively. A typical convergence history 
in terms of the density residual, %, and lift coefficient, Cl, for an NSU2D run is 
shown in Figure 5.2. Approximately 13 Mwords of memory and 60 minutes of 
Cray Y-MP CPU time is required for the solution. 

5.1. Surface Roughness 

In examining the usefulness and practicality of computational methods 
for analysis of true airfoils, the effects of surface roughness of the input 
geometry are addressed. While many common test cases used by code 
developers are smooth, some even defined by a continuous function, 
geometries studied in aircraft analysis are faken from an existing wing. 
The surface description of the TSRV wing comes from splashes, which are 
molded templates taken of the wing, and these splashes were digitally scanned 
to produce the surface coordinates used in the computations. The data tend to 
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include surface irregularities and possible error in the splashes and digitization. 
Preliminary computations on the TSRV geometry revealed high 
frequency oscillations near the airfoil surface, particularly noticeable in the 
surface Cp distribution. The source of the oscillations was determined to be the 
waviness in the surface geometry description by finding a direct correlation with 
the surface curvature. The curvature is defined (as in Ref. 26) as: 


K = 



(5-1) 


This calculation is applied discretely to the surface co-ordinates using central 
difference approximations: 
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Plotting surface curvature, K, along with surface Cp values reveals the 
correlation (Fig. 5.3). 


This noise in the solution not only degrades the neatness of the plots, but 
leads to questions regarding the impact on the turbulence model and artificial 
viscosity, both of which may be influenced by curvature-induced velocity 
gradients at the surface. Therefore, the airfoil geometry was smoothed 
manually by applying a simple spline smoother in a piecewise manner around 
portions of each element. The new geometry was then gridded and 
computations were performed. The smoothed geometry is visually identical to 
the original, but has greatly reduced waviness and the resulting Cp distribution 
is much smoother, as shown in Figure 5.4. Since the Cp of the smoothed 



geometry is essentially the mean of the unsmoothed case, the surface 
roughness appears to cause no adverse effects, such as additional flow 
separation. However, because the Cp plots are cleaner, the smoothed 
geometry is us»^d throughout this study. 

5.2. Cove Shape 

One advantage of Navier-Stokes solvers is the ability to capture flow 
phenomena around complex shapes including sharp corners and concave cut- 
outs such as those created when flaps are deployed. Various alternative 
solution methods require (or are recommended to be run using) streamlined 
surfaces to be used rather than cove cut-outs in order for the solution to 
converge efficiently. One such code is MSES, developed by Mark Drela, which 
is based on a streamline-based formulation of the Euler equations and an 
integral boundary layer mode|27. Potential flow codes also often require similar 
modifications to the surface geom 3try in high-lift configurations. Although the 
difference appears minor, the influence of cove flows on the flow through gaps 
between elements is gaining attention in the literature. Cove fairings which 
deploy with the flaps to reduce the flow separation in the cove regions may 
appear on future aircraft. 

A comparison is made between the true cove shapes for the TSRV in 
approach configuration and the streamlined coves used with MSES in a 
previous study^^. The two geometries and the surface Cp distributions are 
shown in Figure 5.5. The primary effect of the cove fairings is a slightly 
decreased loading on the slat and main element and an increased loading on 
the flap elements due to the altered cove flows. The change in the location of 
the main-element leading-edge pressure peaks is due to the removal of the "lip" 
on the lower surface of the slat in the streamlined case. The flow structure in 
the vicinity of the main element cove is shown for the two cases in Figure 5.6. 



The streamlined cove fairings increase the extent of lower surface attached flow 
and reduce the circulation on that element, thus diminishing the slat effect on 
trailing elements. Therefore, with less damping due to in''isc'd interaction, the 
trailing element loading increases. In addition, the cove fairings reduce the 
wake thickness slightly, causing a reduction of viscous damping on the trailing 
elements. As shown in Figure 5.7, the main-element wake is thinner for the 
streamlined-cove geometry. As a result of the increased fore-flap loading, the 
pressure-recovery demands increase and the onset of separation moves 
forward, from r\ = 0.924 with cut-out coves to r| = 0.913 with cove fairings. The 
increased region of flow separation is indicated by the flattening of surface 
pressure at the fore flap trailing edge (Fig. 5.5), and causes a thicker wake at 
the r\ = 0.950 station, as shown in Figure 5.7. This trend of increased 
separation is a result of this specific case, but generally, if the contour chosen 
for cove fairings approximates the shape of a streamline of the true flow, the 
impact on the computations is small, on the order of other errors evident in the 
solution. 

5.3. Sweep Theory 

The correlation between 2-D computations and flight test data depends 
highly on the applied sweep correction. Figure 5.8 shows computations using 
simple sweep theory, the local sweep correction, and local sweep plus the 
empirical n-factor correction. Experimental data is also shown to reveal errors 
and the MSES computation from Reference 23 (performed using simple sweep) 
is included to illustrate the itnprovement gained by Navier-Stokes solvers. In 
each case, the an required to match the slat pressure distribution was different. 
For MSES, a„ - 1.2®, while for NSU2D using simple sweep theory, a„ = 2.0®, 
with local sweep theory, On = 3.0®, and with the n-factor, Un « 2.5®. As 
expected, the local sweep theory provides closer matching along the main 



element upper and lower surface, but the more accurate correction to the flaps 
reveals additional overprediction of the flap loadings. The n-factor variation 
reduces the flap overprediction drastically, thus making the bes. correlation, but 
perhaps artificially so, since the method is empirical, and without solid 
theoretical justification. 

Several issues remain to be addressed regarding the sweep correction. 
Simple sweep theory accounts purely for sweep, assuming an infinite, 
untapered, swept wing. The local sweep method improves the correction by 
introducing the effects of wing taper. Using either method, the downwash due 
to the vortex system of the 3-D wing is partially accounted for by the a„ 
correction. The downwash effectively reduces the angle of attack of the flow 
past the airfoil section. In addition, spanwise and chordwise variations in 
downwash exist in true flow, and it is possible that these variations explain the 
apparent validation of the n-factor correlation. Gap size and chordwise 
variation in downwash are thought to be tha primary sources of the 
overprediction of flap loadings. 

5.4. Grid Refinement 

In order to have a smooth variation in cell sizes, initial grids were 
generated with high aspect ratio, viscous cells only near the airfoil surfaces. 
Yet, the dissipation of shear layers in the solution was observed to correlate 
with the expanding grid cell sizes (Fig. 5.9). In order to improve resolution of 
the flow, two types of grid refinement were used. First, an automated solution- 
dependant refinement was carried out using Meshad of the NSU2D package. A 
velocity gradient three times the average gradient was the refinement criterion 
used in an effort to improve resolution of viscous gradients. One region which 
is notably resolved in this manner is the shear layer in the cove regions, as 
shown in Figure 5.10. The second metfjod of grid refinement was a whole new 
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generation of grids using Gridpts and UMESH2D. Specifically, longer wake 
r?gions with more nodes in both a streamwise and tangential direction were 
preserved from the initial hyperbolic grids around the elements The additional 
streamwise nodes promote the smooth topology that was previously obtained 
by shortening ihe extent of high aspect ratio cells. A portion of the wake grid in 
the vicinity of the flap elements is shown in Figure 5.11. The sizes of the three 
meshes tested are listed in Table 5.1. 



Baseline Mesh 

Meshad case 

Wake case 

Total Nodes 

72,975 

118,157 

88.799 

Slat Surface 

229 

327 

214 

Main Element Surface 

399 

604 

399 

Fore Flap Surface 

210 

307 

210 

Mid Flap Surface 

363 

454 

319 

Aft Flap Surface 

213 

300 

214 


Table 5.1 - Mesh sizes for five-element grid comparisons. 


The surface Cp distributions and various velocity profiles for solutions on 
the “baseline", "Meshad", and "wake* grids are shown in Figures 5.12 and 5.13, 
respectively. Clearly, the Meshad case is limited by the grid topology of the 
baseline case. If the initial grid is much more coarse than needed, several 
levels of adaptive refinement are necessary since with each level only one node 
is added between existing nodes in the regions of high gradients. Perhaps an 
improved scheme would include logic to add a number of new nodes between 
two existing nodes based on the size of the local gradient. The results for the 
wake grids show much better resolution of viscous layers, particularly apparent 
in the velocity profiles. Mach contours in the slat wake region are shown in 




Figure 5.14. Clearly, the lack of resolution in the baseline case (Fig. 5.9) 
provides very little gradient to cause mesh refinement, so the Meshad case 
does not resolve the slat wake flow. The wake mesh, however, includes 
numerous cells in the wake region, particularly with fine spacing in the direction 
normal to the wake centerline, allowing excellent resolution of the shear layer. 
The slat wake provides a viscous damping effect on the main element, and this 
reduced loading on the main element reduces the upwash on the slat, due to 
the circulation effect as described by A.M.O. Smith. Similarly, the improved 
resolution of the main- and flap-element wakes increases the extent of the wake 
deficit region, causing increased damping on the flaps, which decreases the 
loading. 

in a rigorous analysis of the wake-mesh result, one might note that the 
overall decrease in circulation does decrease the slat loading, so by definition, 
ttn must be increased to match the slat upper surface pressure of the 
experiment. For this reason, the wake-mesh case was computed with Un 
increased from 3.0° to 3.3°. The effect on the pressure distribution was to raise 
the slat- and main-element pressure levels to match those found using the other 
meshes. Yet, the decrease on the flaps remains, since the loading of flap 
elements are not sensitive to small changes in angle of attack. Therefore, the 
wake meshes are still advantageous for resolving viscous-flow phenomena and 
providing a closer correlation with flight. 

One obvious problem with wake grids is that the high-aspect-ratio wake 
cells are concentrated where the panel solution of Gridpts determined the 
wakes to exist. The advantage of an existing grid is reduced at angles of attack 
for which the wakes are not aligned where the grid is well resolved, though the 
larger number of nodes throughout the viscous region does benefit the solution. 
Generation of additional grids for use at the other angles of attack would be 



time consuming since the wake locations are determined at such an early stage 
in the Gridpts grid generation scheme. Perhaps future grid generators could 
place a fan of high aspect ratio cells in the critical regions of the flow, making 
the grid more versatile. 

5.5. Transition Specification 

In complex multi-element flows, laminar regions along the leading edges 
of the elements help to reduce the boundary layer momentum thickness, 
delaying separation and reducing drag while increasing maximum lift. 
Relaminarization of turbulent attachment-line flows has also been observed in 
tests^-3'20. In order for the computations to best predict the flow, transition 
should be modeled. The Spalart-Allmaras turbulence model, like many others, 
allows the specification of laminar regions of the flow along airfoil surfaces. In 
these laminar regions the source term for the eddy viscosity is turned off. 

To assess the effect of specifying transition, the approach configuration 
was analyzed with laminar flow on the slat from ihe upper surface trailing edge 
to the edge of the lip on the lower surface. Figure 5.15 shows the comparison 
between ihe fuliy turbulent result and that with the laminar slat. Although not 
very evident in the Cp distribution. Table 5.2 shows a 5% decrease in drag, as 
well as a 0.7% increase in I'ft due to the laminar flow on the slat, resulting in a 
6% increase in the critical aerodynamic parameter L/D. The velocity profiles in 
Figure 5.16 show the thinner laminar boundary layer on the slat surface which 
decreases the form drag on that element. The laminar velocity profiles also 
have a lower velocity gradient in the immediate vicinity of the slat surface than 
the turbulent case, thus producing less skin friction drag than the fully turbulent 
flow. The thinner slat wake in the laminar case causes less damping of the 
main element pressure, resulting in the slight increase in total lift. 
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Cl 

Cd 

L/D 

Fully Turbulent 

3.509 

0.04878 

72.05 

Laminar Slat 

3.534 

0.04628 

76.36 


Table b.2 - Force coefficients obtained from the surface pressure integration. 

5.6. Comparison With Flight Data 

The preceding sections illustrate several issues faced in the 
computational prediction of viscous flow, using 2-D analysis tools and 
comparing with true 3-D flow. This section shows the best correlation and 
examines the abilities and limitations of the current methods. 

5.6.1. Landing Configuration 

The 40° flap setting used during landing provides the primary geometry 
used in this analysis. The specific case compared here is for a flight angle of 
attack of 8.12°, M« = 0.24, and Re« = 14.67 million. The best correlation 
between 2-D N-S computations and flight data was obtained with extended- 
wake meshes, using the local sweep correction along with the empirical n-factor 
variation, run at Un = 3.0°. As shown in Figure 5.17, the surface pressures are 
well predicted. The computations indicate the location and level of the second 
pressure peak at the main element leading edge due to the discontinuity where 
the slat nests in cruise. Most of the pressure distribution along the main 
element upper and lower surface is predicted accurately. The pressure 
distribution along the flap lower surfaces is accurately predicted, and though the 
flap upper surface pressures are somewhat off, the levels relative to each other 
are good. Although detailed boundary-layer measurements were not conducted 
in Phase II flight tests, tuft flow observations indicated flow separation on the 
fore flap, also characterized by the flattening of the pressure near the fore flap 
trailing edge. Although the computations do not show the flattening of the 




46 


surface pressure distribution as much as experimental data, separation is 
clearly evident in the computed flow field (Fig. 5.18). The usefulness of CFD 
methods to aid in the understanding of flow structures is further illustrated in 
Figure 5.19, which shows total pressure contours over the aft portion of the 
configuration, indicating the growth and interactions of element wakes. 

Currently, CFD methods are best used to indicate qualitative information 
regarding the flow structure, and quantitative data such as variations due to 
changing flow conditions. However, accurate predictions of specific 
performance characteristics, such as lift and drag forces and stall angles are 
still limited. Yet, the comparisons of normal force coefficient, Cn, shown in 
Table 5.3 indicate improvements in this area. 



Experiment 

Simple sweep 

Local sweep, n 

Slat 

0.2240 

0.2218 

0.2320 

Main 

2.2603 

2.5053 

2.3067 

Fore flap 

0.3345 

0.4025 

0.3725 

Mid flap 

0.3515 

0.3696 

0.3502 

Aft flap 

0.0450 

0.0493 

0.0486 

Total 

3.2153 

3.5484 

3.3100 


Table 5.3 - Cn values obtained from the surface pressure integration. 

Several key sources of error can explain the discrepancies between the 
computations and the flight test results. The major source of error is believed to 
be in the geometry, namely in the gaps and deflections. The shapes of each 
ii . dividual element are accurate, matching templates on the actual aircraft, but 
the gaps, overlaps, and angles which significantly affect the interaction of the 
dements were measured more crudely to construct the computational model. 
Furthermore, the geometry used is that measured on the ground, whereas the 




flight data obviously comes from the lifting geometry which changes with 
blowback due to play in the linkages and aeroelastic deformations. 

Along with the g»-ometry problems, the computational model is limited by 
the sweep theory ;'sed to correlate the 2-D results to 3-D dpta. The local sweep 
theory does attempt to account for sweep and taper, and the n-factor 'iorrection 
may account for other 3-D discrepancies which have been observed 
experime. .tally. However, these corrections are relatively simple compared to 
Ihe complexity and detail of the flow sc’ 'ers, and are believed lo be of first order 
accurav.,/. 

Additional error in the computational results may be attributed to known 
numerical errors. Artificial viscosity which provides stability for the 
computations dissipates gradients. Any lack of grid resolution causes additional 
dissipation. Furthermore, turbulence modeling is rapidly changing since current 
methods are still known to fail in regions of separation. The landing 
configuration of the TSRV geometry has four regions of separated flow due to 
the coves and the fore-flap upper surface separation, all contributing to errors in 
circulation as well as the inviscid interactions between the elements. However, 
the excellent agreement between the computations and 2-D wind tunnel data 
shown in Section 3 would suggest that the error in the codes is much less than 
those due to geometry and 3-D flow effects. 

5.6.2. Take-Off Configuration 

The 15° flap setting used for take-off was also tested using wake-type 
grids along with a local sweep correction and n-factor variation. Free-stream 
conditions for the experimental case used here are characterized by = 0.17, 
Re« = 1 1.85 million, and a = 9.5°. The surface pressures are shown in Figure 
5.20. Again, the correlation is very good, with many of the same explanations 
for error. In particular, the slat gap used for the computations is the same as 



‘hat used in Reference 22, but is larger than the gap specified in a recently- 
acquired model of the geometry. This case brings up an additional difficulty in 
the computations. Using the NSU2D package to generate grids which preserve 
cells to resolve wakes requires the structured mesh regions around each 
element, wfiich are increasingly difficult to work with and to obtain a satisfactory 
mesh as the gap between elements decreases. 

5.7. Sensitivity to Changing Flow Conditions 
5.7.1 . Angle of Attack Effects 

The flow past the approach configuration was computed for a range of 
angles of attack, and the corresponding pressure distributions are plotted in 
Figure 5.21. The flight data are presented h Reference 1. The cases of an = 
-2°, 3°, and 8° correlate roughly to flight angles of attack of 2®, 8®, and 13°, 
respectively. The computations were not matched to specific test points 
because of the additional CPU time this would consume in the iteration to find 
the converged solution which best matches the slat pressure distribution. The 
trends are, nevertheless, clearly illustrated using these test cases. The slat and 
main element pressure peaks increase drastically and the slat loading moves 
forward with increased a. At low a, the separated flow in the slat cove 
suppresses the pressure on the main element leading edge lower surface. At 
higher angles of attack, the separated flow is better contained in the slat cove 
region, and the pressure on the main element lower surface is higher. On the 
flaps, for angles of attack in the range of flight data, the variation in loading due 
to a is small, in agreement with the test results. The trend is a small decrease 
in flap loading due to the main element wake thickness growth with increasing 
a. As in the flight test results, the aft flap pressure distribution is independent of 
c.. The On = 13° case represents a higher loading than obtained in flight. The 
main element upper surface develops a thick boundary layer and the 
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consequently thick wake drastically suppresses the flap loadings. At the higher 
loading, the tendency for upper surface separation increases and inaccuracies 
due to the turbulence model become prevalent. 

5.7.2. Reynolds Number Effects 

The sensitivity of the approach configuration was also assessed for 
changes in free-stream Reynolds number. The results are shown in Figure 
5.22. In the range of flight conditions, roughly from Ren = 10- to 16-million 
(corresponding approximately to Re„ = 12- to 20-million), the effect of changing 
Ren is very subtle. Boundary layer and wake thickness change very little in this 
flow regime. On the contrary, the results for Ren = 1 .5 million show a significant 
effect. This case was chosen to approximate a typical wind-tunnel model at 
standard atmospheric conditions. The change in the pressure distribution 
illustrates the problem of extrapolation from wind-tunnel results to predict flight 
performance. At low Ren, viscous effects have a more significant effect due to 
the thicker boundary layers and wakes, which suppress the loading of trailing 
elements. Also, in experiments at different Rep. the mechanisms of transition 
and separation act differently. 

5.8. Gurney Flaps 

The concept of the Gurney flap which originated In racing car spoiler 
design has found practical application on aircraft wings^^. This small flat plate 
placed at the trailing edge of an airfoil, normal to the lower surface, effectively 
shifts the location of the Kutta condition, increasing the circulation and hence 
the lift of the element The increase in lift comes primarily from the effective 
increased camber on the lower surface without adversely disturbing the upper 
surface flow. The profile drag of the airfoil generally increases with the addition 
of the Gurney flap, but often the percentage increase in lift is greater, resulting 
in an increased lift-to-drag ratio and therefore better efficiency and 
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performance. Because of the simple geometry of the Gurney flap, construction 
is simple and weight is low, and deployment witii the flap system is easily 
accomplished. Various theories and experiments describe a flow structure 
downstream of a Gurney flap which has dual recirculation regions^s-^o. 

Gurney flaps of length 0.01c and 0.02c were added to the aft flap in the 
computational geometry. In order to later compare with the results of the 
upcoming Phase III flight tests, the devices were modeled as “U'-shaped 
brackets, including the base required for mounting to the wing lower surface. 
The flow structure in the vicinity of the Gurney flap is shown in Figure 5.23. The 
dual recirculation regions downstream of the flap are predicted as described in 
the literature. No grid refinement was necessary to reveal this phenomenon. 
Some separation e'^ists in the region upstream of the Gurney flap, but the 
primary flow turns smoothly downward with the added camber of the lower 
surface. The Mach contours downstream of the Gurney flap show the shrinking 
of the recirculating wake as the flow moves away from the airfoil (Fig. 5.24). 
The impact of the Gurney flaps on the loading of the system is shown in Figure 
5.25. The circulation of the entire system increases with the addition of the 
trailing*edge device, illustrating the elliptic behavior of the flow and the 
increased circulation of the airfoil system. Along with the increase in lift, the 
Gurney flaps produce a substantial drag rise, primarily in the form of pressure 
drag, as might be expected from the orientation of the Gurney flaps in the 40° 
flap setting. The corresponding lift-to-drag ratio is decreased, as shown in 
Table 5.4. Yet, in the approach configuration, drag can be a favorable quality, 
assisting the deceleration of the aircraft, and high lift is often more important 
than high L/D. 



Cl 

Cd 

L/D 

Baseline Airfoil 

3.509 

0.04878 

72.05 

1®/o Gurney Flap 

3.655 

0.05245 

69.69 

2®/o Gurney Flap 

3.710 

0.05566 

66.65 


Table 5.4 - Force coefficients integrated from the surface pressure 
distributions for the approach configuration (40° flaps). 


The effects of 1% and 2% Gurney flaps on the take-off configuration are 
shown in Figure 5.26. Again, the primary effect of the flaps is an increased 
circulation on the entire system. The corresponding force coefficients are listed 
in Table 5.5. 



Cl 

Cd 

17D 

Baseline Airfoil 

2.208 

0.03615 

61.10 

^% Gurney Flap 

2.299 

0.05419 

42.43 

2 % Gurney Flap 

2.388 

0.06175 

38.67 


Table 5.5 - Force coefficients integrated from the surface pressure 
distributions for the take-off configuration (15® flaps). 
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Figure 5.3 - Correlation between waviness in surface pressure 
computations and airfoil surface curvature, K. 





Figure 5.4 - Effect of smoothing airfoil surface coordinates on computed 

surface pressure distribution. 
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Figure 5.5 - Effect of cove shape modeling on surface pressure. 
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Figure 5.7 - Effect of cove shape modeling on boundary-layer 

velocity profiles. 





Figure 5.3 - Effect of sweep theory on surface pressure distribution. 
(M„ a 0.24, Re, a 14.67 million, a a 8,3®) 




Figure 5.9 - Mach contours overlaid on baseline mesh in vicinity of 
main-element leading edge showing dissipation of slat wake. 

NOTE: This figure was originally in color and the image shown 
here lacks the clarity provided by the color overlay. 
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a. Baseline mesh. 



b. Wake mesh. 


Figure 5.11 > Grid topology in vicinity of fiap elements. 
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Figure 5.13 * Effect of mesh on boundary-layer velocity profiles. 
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a. Adaptively-refined mesh. 



b. Wake mesh. 


Figure 5.14 - Mach contours overlaid on computational mesh in virinity of 
main-element leading edge showing resolution of slat wake. 

NOTE: This figure was originally in color and the image shown 
here lacks the clarity provided by the color overlay. 
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Figure 5.16 • Effect of laminar slat ilow on boundary-layer velocity profiles. 





Figure 5.17 - Best correlation with flight results for approach configuration 
(40° flaps), using wake meshes, local sweep correction -i- n-factor. 

(M„ = 0.24, Re^ = 14.67 million, a = 8.3°) 



Figure 5.18 - Velocity vectors showing separated-flow region on 
fore-flap trailing-edge upper surface. 
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6. Conclusions 

The impact o* high-lift performance on the overall performance of a 
subsonic transport aircraft justifies comprehensive studies of the flow past multi- 
element airfoil systems. Numerous wind-tunnel- and flight-test programs have 
provided tremendous insight into the high-lift flow regime, bsit only with the 
support of a tremendous budget. In an effort to reduce the expense of 
analyzing numerous configurations and the effects of changes in freestream 
conditions and geometry, recently developed Computational Fluid Dynamics 
(CFD) methods are being applied to multi-element airfoil systems. As part of a 
multi-phase high-lift flight research program, computational studies are being 
conducted on the NASA Langley Transport Systems Research Vehicle (TSRV). 
In this study, the capabiliiies of two production codes which solve the 2-D 
Navie -Stokes equations on unstructured grids are assessed. Computational 
results are compared to wind-tunnel and flight-test data. The effects of surface 
roughness, cove shape modeling, sweep theory, grid refinement, and transition 
specification are addressed. In addition, the trends associated with charging 
freestream angle of attack and Reynolds number are shown. Finally, the 
aerodynamic effects of Gurney flaps are predicted. Several primary 
conclusions resulting from this study are summarized here; 

• CFD tools accurately predict the 2-D viscous flow past multi-element wing 
sections tested in wind tunnels in the absence of massive upper-surface 
separation. 

• Application to 3-D flight-test studies introduces additional complications in 
terms of geometry inaccuracies and sweep theory to correct 2-D 
computations for the effects of sweep, taper, and finite wing effects. 
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• Surface roughness in the airfoil geometry used for computations causes 
waviness in the surface pressure distribution, but does not c ^stically affect 
the solution in ter ns of the average levels, or the onset of separation. 

• The common practice of modeling coves with smooth fairings for 
computations causes a small decrease ’n circulation on the altered 
elements, and slightly thinner wakes, and therefore a slightly increesed 
loading on trailing elements. However, the overall effect on the flow is small. 

• Sweep theory has a first-order influence on the comparison with 3-D flight 
data. An improvement over simple sweep theory accounts for taper effects 
by using local sweep angles, which vary along the chord, rather than a 
nominal angle. An empirical "n-factor" alteration to the sweep correction 
applied to Cp improves the comparison with flight, but without solid 
theoretical justification. 

• Grid refinement improves resolution of viscous-flow phenomena, such as 
wakes and confluent boundary layers, but with only a small effect on the 
surface pressure distribution. 

• Generating grids with additional cells in wake regions, and in particular, high 
aspect-ratio cells with fine spacing normal to streamlines is a more efficient 
means of resolving viscous flows than using adaptively-refined meshes, 
which require several levels of refinement and introduce more new nodes. 

• Specifying a region of laminar flow using the Spalart-Allmaras turbulence 
model causes insignificant change to the flow prediction. This result may 
seem to show a weakness of the model, but it may also be interpreted as an 
indication that the fully-turbulent solution cotrectly applies small levels of 
eddy viscosity in regions where the true flow is laminar. 

• Using the locai sweep correction, along with the n-faclor modification, the N- 
S computations using wake meshes correlates well with flight-test results. 



including the prediction of sepa-ated-flow regions. Discrepancies are 
attributed to sweep theory and inaccuracies in the model geometry. 

• MSES provides an excellent prediction of surface pressures using much 
less memory and CPU time than N-S solvers, but is more limited in terms of 
geometry and flow conditions, and provides less detail of the flow structure. 

• The computed trends of changing angle of attack agree with flight results. 
Slat and main-element loading increases with increasing a, while flap 
loading decreases slightly due to the slat effect as well as viscous damping 
resulting from the thicker wakes. However, due to known inaccuracies in 
the turbulence modeling, the prediction of maximum lift was not attempted. 

• The computations show negligible sensitivity of the TSRV wing section to 
changes in Reynolds number in the range of flight. A prediction using a low 
Reynolds number, comparable to that of an atmospheric wind tunnel, shows 
a significant decrease in circulation and loading. 

• The addition of 1% and 2% Gurney flaps to the aft flap produce the 
expected flow structure, as well as lift increments, but a decrease in L7D for 
this configuration. 

The results of this study illustrate the capability of current computational 
methods to accurately predict the flow past multi-element airfoil systems. 
Furthermore, 2-D solvers, which are manageable with current computer 
resources, can provide an accurate representation of 3-D flow, when used in 
conjunction with sweep theory. The sweep corrections, applied here in the form 
of pre- and post-processing operations, attempt to account for the effects of 
wing sweep, taper, and downwash from the tip vortices. Until 3-D viscous-flow 
solutions become practical, further improvements to the intermediate step of 
quasi-3-D solution methods are necessary. A potential improvement would be 
to add the third dimension to the equations, but still solve the equations on a 



73 


2-D grid, with gradients in the direction normal to the 2-D cut turned off. This 
would eliminate the need for pre- and post-processing operations and would 
allow the development of a 3-D boundary layer. Another area in need of 
attention, particularly in the correlation with flight test data, is the modeling of 
the geometry. Multi-element airfoil flows are highly sensitive to changes in 
incidence angles and gaps. The geometry measured on the ground may be 
quite different from that in flight under loading, so the computational model 
should, if possible, be made from the loaded system. For this reason, 
upcoming flight tests of the TSRV will include OPTOTRAK instrumentation to 
measure movement occurring in flight. Despite current limitations and known 
errors, CFD methods can provide tremendous insight of multi-element high-lift 
flows. 



Appendix 1: Upwinding of Convection Terms 


CFD methods require smoothing or diffusion for stability and 
consequently to converge on a solution. This diffusion is commonly aoded 
explicitly into the equations being solved. An alternative method of achieving 
the same result is upwinding of the convective terms in the discrete formulation. 
While the upwind method applies to multi-dimensicnal flows solved using finite- 
element and finite-volume methods, the analysis of this method is most clearly 
shown for one dimension using a finite-difference approach. The model 
equation used nore is the simple linear convection equation, known as the wave 
equation: 


du _ du 

Tt~~^Tx 


(AM) 


If artificial viscosity was explicitly applied, .he equation solved would take 
the form, 

(A1-2) 


du du d^u 

~dt~~^lx^^'dj? 


where a is the coefficient of artificial viscosity, a« 1. Note that a must be a 
function of Ax in order for Eq. (A1-2) to be consistent with Eq. (A1-1) as Ax-»0. 

Using upwind differences, the convective terms of Eq. (A1-1) are 
discretized as: 



-c- 


M. -Mi. 


-c 


Ax 
Ax 


,c>0 
, c<U 


(A1-3) 


Applying truncation error analysis to the differenced terms In Eq. (A1-3), we 
expand the terms that are not at time level, n, and at space index, I, using 
Taylor series expansions. For the case when c> 0, we have: 
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Wi_, =«, 


A.r‘ 


2 (?jc' 


-O(Ax^). 


(A1-4) 


In a complete truncation error analysis, the discretization of the time-dependent 
derivative would also be evaluated, using Taylor series expansions of those 
terms v'hich are not at time level, n. However, the purpose of this discussion is 
to show the relationship between different discrete forms of the space 
derivative, so it is assumed that the same time discre.ization is used in both 
cases, and an analysis of the discrete time terms is not shown here. 

Substituting Eq. (A1-4) into Eq. (A1-3), for c> 0 yields: 


-c 


it. 

dx 


/(" - u" + /^Jr— 

Ax^ d^u 

+ 0{Ax^) 

i 

/J, u. 

, Td? 

A.X 


(A1-5) 


which reduces to: 






(A1-6) 


As might be expected from this one-sided difference, the truncation error of the 
discrete formulation is first order. Ignoring the terms of second order and 
higher, the numerical error, e in the upwind difference method shown here is: 


e = 




(A1-7) 


Thus, we see the artificial viscosity, a is evident in the solution, as if it were 
jxplicitly applied as in Eq. (A1-2). For the case when c < 0, the leading error 
term is: 


e 



(A1-8) 


We confirm that upwinding is a special case of adding artificial viscosity 
to a central difference scheme by observing the discrete form of the right hand 
side of Eq. (A1-2): 
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Substituting a~c — yielcis: 
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2 J Ax^ 




Ax 


(AMO) 


which is identical to tne upwind difference in Eq. (A1-3) for o > 0. Similarly, 

Ay 

substituting a = ''“'to (A1-9) yields the upwind difference for c < 0. 



Appendix 2: Local Sweep Angles 


For a tapered wing, the angle of sweep, A, varies along the chord. 
Although one might assume a linear variation in A from Ale Ate. if is actually 
tanA that varies linearly along the chord, as shown hore. This geometric 
derivation is for a single-element wing. In order to find local values of A on a 
multi-element wing. Ale and Aje must be known for each elemetit. 

Figure A2.1 shows a generic swept, tapered wing. The root chord length 
is Cr, so for a taper ratio, X, the tip chord length is Ct = XCf. The wing semispan 
is b, and a is the distance in the x-direction to the wing tip leading edge due »o 
the leading-edge sweep. From these dimensions, we find relations for the 
leading-edge and trailing-edge sweep angles; 


tanAuj = 


a 


(A2-1) 


tanA 


re 


_('! + Ac,.) - c, + CfjX - 1) 


b b 

Solving Eq. (A2-2) for Cr and substituting Eq. (A2-1) for a yields: 


(A2-2) 


^tan Are -fctan A,j 

(T^li) 


(A2-3) 


The chordwise location on the airfoil, non-dimensionalized between the leading 
edge and the trailing edge, is represented by For an arbitrary chord ..lation, 
the local sweep angle, A^, is given by: 


^ b b 

Substituting Eq. (A2-1) for a and Eq. (A2-3) for cr yields: 

tan A^ - 


(A2-4) 


b 


(A2-S) 



and simplifying, we find: 


ian = (1 - (^)tan A,jr + ^lan (A2-6) 

Thus, the tangent of the local sweep varies linearly between the tangent of the 
leading-edge sweep and the tangent of the trail'ng-edge sweep. 
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